Surface diffusion and low vibrational motion with interacting adsorbates: 
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Here, an approach in terms of shot noise is proposed to study and characterize surface diffusion and 
low vibrational motion when having interacting adsorbates on surfaces. In what we call statistical 
limit, that is, at long times and high number of collisions, one expects that diffusing particles display 
an essential Markovian behavior. Accordingly, the action of the pairwise potentials accounting for 
particle-particle collisions is equivalent to considering a shot noise acting on a single particle. We 
call this approach the interacting single adsorbate approximation, which gathers three important 
advantages: (i) the dynamics underlying surface diffusion and low vibrational motion can be easily 
understood in terms of relatively simple stochastic processes; (ii) from our model, appropriate (and 
well justified) working formulas are easily obtained, which explain the results arising from more 
complicated (but commonly used) molecular dynamics simulations within the Langevin formulation; 
and (iii), at the same time, it is less demanding computationally than the latter type of calculations. 
In order to illustrate the application of this model, numerical results are presented. Specially, our 
model reproduces the experimental observation regarding the broadening of the quasielastic peak 
ruling surface diffusion. 

PACS numbers: 05.10.Gg, 05.40.-a, 68.43.-h 



I. INTRODUCTION 

Diffusion of adsorbates (e.g., atoms or small molecules) 
on metal surfaces is one of the most fundamental pro- 
cesses in surface science for which there is a wealth of 
available experimental data. These data are obtained 
by means of either standard time-of-flight techniques 

0, H n 

or the more novel spin-echo measurements [J]. 
In both cases the theoretical approach followed to under- 
stand and explain the experimental results is basically 
the same. At low coverages @ adsorbate-adsorbate inter- 
actions can be neglected; diffusion {self- diffusion) is well 
characterized by only studying the dynamics of a single 
adsorbate interacting with the surface. This approach 
is known as the single adsorbate approximation, where 
adsorbate dynamics are described by means of the (clas- 
sical) standard Langevin equation, i.e., the diffusionpro- 
cess is considered as a Brownian-like motion [l|, H, Hf, 0] • 
On the other hand, as the coverage increases, adsorbate- 
adsorbate interactions can no longer be neglected. In 
these cases, molecular dynamics techniques within the 
Langevin framework (MDL) are commonly used to study 
the problem to deal with the surface thermal vibrations. 
In this type of calculations the lateral dipole-dipolc inter- 
action between Na atoms as a function of the coverage is 
accounted for by pairwise potential functions. Previous 
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works [3, [8| using this approach have failed when trying 
to reproduce the experimental data quantitatively. How- 
ever, recently it has been shown [9( that a good agreement 
with the experiment can be achieved by considering the 
adsorbate motion in three dimensions. 

Introducing pairwise potential functions results in a re- 
alistic description of the adsorbate dynamics. However, 
there is not a simple manner to handle the resulting cal- 
culations by means of a theoretical model, and one can 
only proceed by using some suggested fitting functions. 
Moreover, employing MDL simulations always result in a 
relatively high computational cost due to the time spent 
by the codes in the evaluation of the forces among par- 
ticles. This problem becomes worse when working with 
long-range interactions, since a priori they imply that one 
should consider a relatively large number of particles in 
order to obtain a good simulation. In such cases, one can 
also make use of truncations, as is commonly done when 
considering short-range interactions, but they might lead 
to serious inaccuracies (Toj |. Some alternative techniques 
have been proposed to overcome this drawback, though 
they are more expensive computationally. 

In order to avoid such inconvenient (interpretative and 
computational) one can try to translate the behaviors ob- 
served in the experiment into a simple, realistic stochastic 
model. To achieve this goal, first note that when a few 
adsorbates are present on a surface, the single adsorbate 
approximation perfectly describes their dynamics [ill ]. 
For practical purposes, this means that long-range inter- 
actions can be assumed as negligible; otherwise, impor- 
tant effects would appear even when having adsorbates 
located at very far distances from each other. Hence, in 
principle, when having a higher coverage, one can assume 
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that one adsorbate will basically feel the action of other 
adsorbates located in a certain neighborhood (as is cor- 
roborated by the calculations presented in Ref . , where 
only about 400 atoms were used in the simulations). In 
other words, the importance of the long-range tail of the 
particle-particle interaction becomes only relative. Now, 
in order to obtain a good simulation of a diffusion pro- 
cess, one has to consider very long times in comparison 
to the time scales associated to the friction caused by the 
surface or to the typical vibrational frequencies observed 
when the adsorbates keep moving inside a surface well. 
For example, in the Na/Cu(001) system the time scales 
associated to typical vibrational frequencies and frictions 
are of the order of a few picoseconds, while the propaga- 
tion time employed in a typical MDL simulation is of the 
order of several nanoseconds [9( . This means that there 
will be a considerably large number of collisions during 
the time elapsed in the propagation, and therefore that 
at some point the past history of the adsorbate will be ir- 
relevant. That is, adsorbates undergo so many collisions 
after some time that they lose any trace of the type of 
interaction among them when they are considered statis- 
tically. This memory loss is a signature of a Markovian 
dynamical regime, where adsorbates have reached what 
we call the statistical limit. Otherwise, for time scales 
relatively short, the interaction is not Markovian and it 
is very important to take into account the effects of the 
interactions on the particle and its dynamics (memory 
effects). 

Based on the previous physical considerations, here we 
propose a numerical and analytical scheme to deal with 
the diffusion problem of interacting adsorbates as a func- 
tion of the coverage. In particular, our scheme is framed 
within a Langevin formulation, and inspired by the the- 
ory of spectral-line collisional broadening developed by 
Van Vleck and Weisskopf [l2| and the elementary kinetic 
theory of gases 13]. We call this approach the interact- 
ing single adsorbate approximation. The diffusion of a 
single adsorbate is modeled by a series of random pulses 
within a Markovian regime (i.e., pulses of relatively short 
duration in comparison with the system relaxation) sim- 
ulating collisions between adsorbates. In particular, we 
describe these adsorbate-adsorbate collisions by means of 
a white shot noise [14j as a limiting case of colored shot 
noise. In this way, a typical MDL problem involving 
N adsorbates is substituted by the dynamics of a sin- 
gle adsorbate, where the action of the remaining N — 1 
adparticles is replaced by a random force given by the 
shot noise. This type of noise has also been applied, for 
example, to study thermal ratchets [l5j and mean first 
passage times [16| . A different type of collisional model 



has been recently used to obtain jump distributions 17 1. 
A general account on colored noise can be found in [18|. 

Note that since our model faces the problem of adsor- 
bate diffusion from a stochastic perspective, it is alter- 
native to those other models that imply the use of MDL 
simulations. There are no incompatibilities between both 
approaches, stochastic and deterministic ("determinis- 



tic" in the sense of how adsorbate-adsorbate interactions 
are handled), other than their respective ranges of valid- 
ity. On the contrary, they complement each other, since 
the type of results rendered by our model could be used 
to better understand the behaviors and trends displayed 
by the fitting functions that appear in the latter. 

The main goal in this work is to show that a lot of in- 
formation about diffusion can be obtained by means of a 
simple stochastic model based on a shot noise. Thus, no 
substrate friction will be considered in our calculations. 
The combined action of two noises, Gaussian white noise 
to simulate the surface temperature and the white shot 
noise with a nonseparable interaction potential are tack- 
led in a separate publication [l!| (Hanggi and co-workers 
[l5| have also used this combined scheme in the context 
of Brownian motors). In Sec. |TT] we thus give a detailed 
account of the model that we propose here to simulate 
the interacting particle dynamics. Moreover, a discussion 
on the validity of the model is also provided in the light 
of the fluctuation- dissipation theorem. A brief account 
on the broadening of the quasielastic peak with cover- 
age is provided in Sec. IIIII Numerical results for both 
flat and corrugated periodic surfaces (separable interac- 
tion potentials) are shown in Sec. HVI as a function of the 
coverage, temperature, and AA". Finally, the main con- 
clusions extracted from this work as well as a discussion 
on the range of validity of our model are given in Sec. [V] 



II. INTERACTING ADSORBATES AND SHOT 
NOISE 

Particle motion under two-dimensional (2D) separable 
interaction potentials can be simplified as two indepen- 
dent ID motions provided the correlations between the 
two components of the noise source along each direction 
can be neglected. Then, the motion of an adsorbate sub- 
jected to the action of a bath consisting of another ad- 
sorbates on a static 2D (separable) surface potential can 
be well described by a generalized Langevin equation 

x(t) = - f ~f(t-t')x(t')dt' + F[x{t)}+5R{t), (1) 
Jo 

where x represents any of the two adsorbate degrees of 
freedom. In Eq. ([I]), 7(i) is the bath memory function; 
F = — W is the deterministic force, V(x) = V(x + a) 
being a deterministic, phenomenological adiabatic po- 
tential accounting for the adsorbate-surface interaction 
at T = 0, with period a along the x direction; and 5R(t) 
is the stochastic force fluctuation, defined as 



with 



8R(t) = R(t)-((R)), 



((R^^Y.PkWW)^- 

K 



(2) 



(3) 



In this last expression, the double average bracket in- 
dicates the average over both number of collisions (K) 
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according to a certain distribution, Pk, and time (T); 
the subscript T denotes the average along the time in- 
terval T. Note that both F and 5R are forces per mass 
unit. 

The random force R(t) has the functional form of a 
shot noise, 



K 



R(t) = J2b k (t-t k ), 



(4) 



fe=l 



since it describes the collisions among the adsorbates. 
The information about the shape and effective duration 
of the fcth adsorbate-adsorbate collision at tk is thus pro- 
vided by bk(t — tk). According to the definition of a shot 
noise, the probability to observe K collisions after a time 
T follows a Poisson distribution [3] , given by 



Pk{T) 



(AT) 
Kl 



K 



-XT 



(5) 



where A is the average number of collisions per time 
unit. Assuming sudden adsorbate-adsorbate collisions 
(i.e., strong but elastic collisions) and that after-collision 
effects relax exponentially at a constant rate A', the 
pulses in Eq. (fj| can be modeled as 



6 fc (t - t fc ) =c fc A'e- y( *- tfc) , 



(6) 



with t — tk > and Ck giving the intensity of the collision 
impact. Within a realistic model, collisions take place 
randomly at different orientations and energies. Hence it 
is reasonable to assume that the Ck coefficients follow an 
exponential law, 



/(cfc) 



f 



c k > 0, 



(7) 



where the value of a will be determined later on. It can 
also be easily shown that this kind of distribution renders 



the same results as considering the pulse intensity having 
the same value, C, for any collision. In such a case, 
wherever (ck) and (c 2 ,) appear they have to be replaced by 
C and C 2 , respectively. Finally, we would like to mention 
that numerical tests using a different shape function (in 
particular, a Gaussian function) gives the same type of 
results. 



Independently of their intensity, it is apparent from 
Eq. ([6]) that any pulse decays at the same rate A'. This 
rate defines the decay time scale for collision events, r c = 
1/A'. If r c is relatively small (i.e., collision effects relax 
relatively fast), the memory function in Eq. |T]) will be 
local in time. This could be the case, for instance, of 
relatively dilute systems, where the average time between 
consecutive collisions is long enough in comparison to the 
energy transfer process that occurs during the collision. 
Then j(t — t') ~ 7<5(i — t') and the upper time limit can 
be extended to infinity. In doing so we obtain a standard 
Langevin equation 



x(t) = -\x(t) + F[x(t)} + 6R(t). (8) 
Here the friction coefficient 7 measures the number of 
collisions per time unit, as A in Eq. ([5]). Hence, from 
now on, we replace 7 by A. The collisional friction co- 
efficient introduces a time scale r r = I/A, which can be 
interpreted as the (average) time between two successive 
collisions. Thus, although each individual collision lasts 
for a time scale given by r c , its effects over the system 
take place in a time of the order of t t in getting damped. 

In order to obtain some relevant information about the 
adsorbate diffusion process, it is important first to ana- 
lyze the solutions of Eq. JSJ) and then to derive some 
average magnitudes of interest. The former are straight- 
forwardly obtained by formal integration, this yielding 



J 



v(t) = v e- xt + / e- A(t - i ' ) F[x(0]^'+ / e- xit - t ' ) SR(t , )dt l , 
Jo Jo 



x(t) - x + ^(I-e- A *) + i 



t r 



1 - e 



-A(i-t') 



F[x{t')} dt' 



1 

A Jo 



t r 



I - e 



-A(i-t') 



(9a) 

SR(t')dt', (9b) 



r 



where vq = v(0) and xq = x(0). As can be seen, for SR 
I). Eqs. ([9|) are the formal solution of purely (dissipative) 
deterministic equations of motion. Hence, without loss 
of generality, they can be more conveniently expressed as 



v(t) 
x(t) 



v d (t) + v s (t), 
Xd(t) + x s (t), 



(10a) 
(10b) 



and s those other depending on the stochastic force. Note 
that when SR ^ the deterministic part will also present 
some stochastic features due to the evaluation of F(x) 
along x{t). 



where d embraces the deterministic terms of the solutions 
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Since (6R{t)) = 0, 



(v(t)) = v d (t), 

(x(t)) = x d (t), 
(x 2 (t)) = xl(t) + (x 2 s (t)), 



(11a) 
(lib) 
(11c) 
(lid) 



where the barred magnitudes indicate the respective av- 
erages of the deterministic part of the solution, and 



(«?(*)} = e 



-2Xt 



dt' e 2At ' 



t-t' 



(x 2 s (t)) = ^ y o df 



l_ e -A(t-t') 



1 _ g-A(t-i'-r) 



r G(T) dr, (12a) 



G(t) dr. (12b) 



In these expressions, Q (r) is the time correlation function 
of the stochastic force, 

Q(t) = ((SR{t)SR(t'))) = ((SR(t)SR(t + r))) (13) 

[the double bracket is defined as in Eq. ((Sj)]. As seen 
in Eq. (fT2"| . Q{t) is the same for both averages despite 
the time-dependent prefactors being different. Averaging 
over time implies that Eq. (|13p is independent of time; 
average values and time-correlation functions will not de- 
pend specifically on the origin of time, but on the differ- 
ence r between two times considered. This condition 
defines the process as stationary. 



A general expression for Q (r) is readily obtained after 
straightforward algebraic manipulations to be 

G(t) = ^ E W)K [ T {b(t - t')b(t + T- t')) c dt', 
1 K Jo 

(14) 



where 



/(c) dc 



(15) 



is the average over the pulse intensity, with g(c) given by 
Eq. ©. Taking into account that Y iK P k(' t ') k = ^ 
and introducing the change of variable C = t—t', Eq. |T] 
can be approximated by 



G(r) = X (b(0b(C + r)) c dC, 



(16) 



which is a general expression independent of the pulse 
shape. This approximation relies on the hypothesis that 
6(C) ~ outside a narrow time interval < ( < A, 
with A a few times larger than r c , but of the same order 
of magnitude. In particular, substituting Eq. © into 
Eq. (JTSD leads to 



G(r) 



XX' 



-A It 



(17) 



By changing the order of the integration variables, in- 
tegrating over t', and taking advantage of the property 
Q{—t) — Q(t), Eqs. (TT2")) can be expressed as 



1-e 



2_ f X' 

a 1 \ 2(X' — A) 
2 (t 2A'-A 
^2 \a ~ (A' - A)A 2 



-i\t 



X'X 



X' 2 - X 2 



1 



(A' - A) A' 



1 



1-e 



-A't 



-At 



1 - e 
A' 



-(A' + A)i 



2(A' - A)A 2 
1 



A' 2 - A 2 



1 - e 

(A' + A)t 



-2At 



(18a) 



(18b) 



In order to determine the value of a, we assume that 
A' » A, this rendering 



(v 2 (t)) 



a 2 V 



-2At 



2 A 2 



2A< + 1 - f 2 - e 



,-At 



(19a) 
(19b) 



Moreover, when t — * oo, we also assume that the equipar- 
tition theorem holds, and therefore 



1 



m(v 2 (oo)) 



1 



kuT. 



(20) 



2 x y " 2 
Taking into account that v^t) = v 2 e~ 2Xt and that the 



time-dependent term in Eq. (|19a[) vanish asymptotically, 
we obtain a — y/m/ksT. On the other hand, if we 
consider that the system is initially thermalized (i.e., it 
follows a Maxwell-Boltzmann distribution in velocities) 
and has a uniform probability distribution in positions 
around x = 0, then vq = 0, v 2 = hsT/m, and xq = 0. 
Thus, for A' 3> A (i.e., in the Poissonian white noise 
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limit), Eqs. (jlip become 
(«(*)) = o, 



fceT 



<^ 2 W) = 

7 

(x(t)) = 0, 



(x 2 (t)) = xl + 



k B T 
mX 2 



2At+l- (2-e 



(21a) 
(21b) 
(21c) 

(21d) 



These equations constitute a limit. Therefore, for values 
of the parameters out of the range of the approximation, 
deviations are expected. As will be seen, the behavior of 
the numerical results presented in Sec. IIVI fit fairly well 
the trends given by these equations, though A' is not 
much larger than A purposedly. For example, for high 
values of A, the equilibrium thermal velocity is given by 
2//a 2 , with / = A'/2(A' + A), instead of 1/a 2 [note that 
if A' = rA, then / = r/2(l + r), which approaches 1/2 
when r — > oo] . 

After using the approximation A' ^> A, Eqs. (|2"Tj) co- 
incide with those obtained for a Brownian motion. This 
is because within this approximation the shot noise be- 
haves as a (Poissonian) white noise, displaying a behavior 
analogous to that of a Gaussian white noise (Brownian 
motion) whenever the number of collisions per time unit 
(A) is very high and/or the total propagation time (T) 
considered is sufficiently long. Note that in the limit 
where the number of collisions goes to infinity, the Pois- 
sonian distribution approaches a Gaussian one because 
of the central limit theorem. 

As happens for pure Brownian motion (V = 0), two 
regimes are clearly distinguishable from Eqs. (|2ip when 
comparing A and t. For Xt <C 1, collision events are 
rare and the adparticle shows an almost free motion with 
relatively long mean free paths. This is the ballistic or 
free-diffusion regime, characterized by 



(x 2 (t)) 



k R T 



■t l 



(22) 



On the other hand, for Xt 3> 1, there is no free diffu- 
sion since the effects of the stochastic force (collisions) 
are dominant. This is the diffusive regime, where mean 
square displacements are linear with time, i.e., they fol- 
low Einstein's law, 



(x 2 (t)) 



2k B T 
mX 



t = 2Dt. 



(23) 



In analogy to Brownian motion, from this last expression 
we also note for systems driven by a shot noise that (1) 
lowering the density of the particle gas (or, equivalently, 
A) leads to a faster diffusion (the diffusion coefficient D 
increases), and (2) the latter becomes more active when 
the gas is heated. 

Apart from the averages given above, it is also mean- 
ingful to compute the velocity autocorrelation function, 



C(t) = (v(0)v(t))= lim i 



T 



v(t) v(t + t) dt, (24) 



where the correlation time, which provides information 
about the line shape broadening (see Sec. IIIip . is 



1 f°° 

W)L c[T)dT - 



(25) 



The velocity autocorrelation function can also be ex- 
pressed as 

C(r) = {v(t)v(t + T)) 

ft ft+T-t' 

= e -2\t-\r dt , e 2\t> I e Xs g(s)ds. (26) 

JO J-t' 



Introducing Eq. (jT7|) into the right hand side of the 
second equality of Eq. (|2"6")l . and in the limit of large t, 
we reach 



C(t) 



k B T X' 2 X I e- XT e~ x ' T 



A' 2 -A 2 \ A A' 



which in the limit A' 3> A becomes 
C ( T ) = —zr 



At 



. (27) 



(28) 



This expression has been derived starting from a col- 
ored shot noise. However, observe that it is identical to 
the velocity autocorrelation function corresponding to an 
Ornstein-Uhlenbeck process, the only stationary Gaus- 
sian diffusion process. This is because of the approxima- 
tion used here: the number of collisions per time unit is 
very large. In this case, the central limit theorem applies 
and, according to Doob 's theorem, the corresponding cor- 
relation function will display an exponential decay. 

For V ^ 0, Eqs. (|21~j) are expected to show some devi- 
ations at short time scales because of the role played by 
the terms corresponding to the deterministic force. This 
force leads to the presence of a confining potential, which 
makes the particles to move within a bound space region 
for a time. Hence the system will be more localized than 
in the free-potential case (V = 0) and there will be much 
less diffusion. Moreover, some in-phase correlations will 
also be apparent due to the particles that do not have 
energy enough to overcome the potential energy barrier. 
This gives rise to a nondiffusive type of motion. That 
is, a low frequency vibrational motion exerted by the so- 
called frustrated translational mode or T mode. 

The velocity autocorrelation functions given by 
Eqs. (|27|) and (|28|) correspond to a dynamics ruled by 
completely flat surface. Physically, this would be the 
case of a low corrugated surface, for which the static force 
is negligible in both directions. On the contrary, when 
the surface corrugation is important, effects mediated by 
the T mode are expected to manifest in the correlation 
function. An interesting example to examine is that of a 
harmonic oscillator, which can model at a first order of 
approximation the oscillating behavior associated to the 
T mode. In this case d, [H 



C(r) 



k B T 



X 

2uj 



— e- XT/2 ( cos lot - — sin lot I , (29) 
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where 



oj = \ ui. 



A 2 



(30) 



and Wo is the harmonic frequency; for a nonharmonic 
potential, Uq represents the corresponding approximate 
harmonic frequency. Equation (|29p can also be written 
as 



C(r) 



kBT ^e-^ cos (u,T + 5), (31) 



with 



S = (tan)- 1 (^-Y 



(32) 



Observe that Eq. ([28]) is easily recovered in the limit 
ujq -> from either Eq. ((29j) or Eq. (|3l|). 

Finally, a brief discussion on the validity of Eq. |(5J) is 
worthy. As said above, this equation can only be used rig- 
orously when having a white noise, i.e., when Q{t) ~ <5(r) 
and the memory function can be substituted by a con- 
stant. Otherwise, with colored noise, one has to use the 
generalized Langevin equation ([T]) |18]. According to the 
fluctuation-dissipation theorem [2l| . the friction in the 
Langevin equation is related to the fluctuations of the 
random force. This is formally expressed by the rela- 
tionship between the frequency spectrum of the memory 
function and the random force correlation function, 



7 (w) 



k B T 



G{r)e 



'dr. 



Introducing Eq. (jXTj) into Eq. (J33J) yields 

A' 



7M = A T7 



A' + iuj ' 



whose real part is 



Re[ 7 (o;)] = i[7H+7»] = A^^ 



(33) 



(34) 



(35) 



Two limits are interesting in this expression: A' <C 
and A' 3> ui. The first limit involves very short time 
scales (smaller than r c ), where the particularities of 
the adsorbate-adsorbate interaction potential have to be 
taken into account (for example, pairwise interaction), 
and therefore a MDL prescription has to be followed. In 
this case, Eq. (|34|) can be written as 



A' 2 

7(w) w A — j , 



(36) 



and hence the standard Langevin equation can no longer 
be used. Moreover, as a consequence, since the frequency 
depends on the friction the relaxation time scale r r is not 
well defined. Conversely, in the second case, the collision 



time scale rules the system dynamics, since it establishes 
a (frequency) cutoff. This leads to 



7 (w) « A 1 



A 72 



(37) 



which can be written as 7(0;) ~ A [notice that A = 7(0)] 
whenever A <C u> <C uj c — r c _1 - As can be seen, in this 
limit (which can also be written as r c <C r r ), the use of a 
standard Langevin equation is well justified. This limit 
holds for strong but localized (or instantaneous) collisions 
(as assumed here) as well as for weak but continuous 
interactions (Brownian motion). 



III. ELEMENTS OF SURFACE DIFFUSION 

In diffusion experiments carried out by means of 
quasielastic helium atom surface scattering (QHAS), one 
usually measures the differential reflection coefficient. In 
analogy to liquids [22j this magnitude is given as 



= n d TS(AK,to) 

= n d T ff G{R,t)e^ AK - R -^ dRdt. 



d 2 K{AK,uj) 
dQdiu 



(38) 

This expression is the probability that the probe (He) 
atoms scattered from the diffusing collective (chattered 
on the surface) reach a certain solid angle Q with an en- 
ergy exchange huj = Ef — Ei and wave vector transfer 
parallel to the surface AK = K/ — K;. In Eq. (f38|) . rid 
is the (diffusing) surface concentration of adparticles; T 
is the atomic form factor, which depends on the inter- 
action potential between the probe atoms in the beam 
and the adparticles on the surface; and S(AK,w) is the 
dynamic structure factor or scattering law, which pro- 
vides a complete information about the dynamics and 
structure of the adsorbates through particle distribution 
functions. The dynamic structure factor is therefore the 
observable magnitude that we are interested in here. Ex- 
perimental information about long distance correlations 
is obtained from it when using small values of AK, while 
information on long time correlations is available at small 
energy transfers, fku. On the other hand, a standard 
procedure employed to obtain the adiabatic adsorption 
potential V^R) mediating the adsorbate-substrate inter- 
action consists in starting with a model potential that 
contains some adjustable parameters. Then, Markovian- 
Langevin equations are solved for different friction coef- 
ficients with Gaussian white noise to reproduce the ex- 
perimental QHAS measurements [2[. 

When dealing with interacting particles, particle distri- 
bution functions are described by means of the so-called 
van Hove or time-dependent pair correlation function 
G(R, t) Given a particle at the origin at some ar- 

bitrary initial time, G(R, t) represents the average prob- 
ability for finding the same or another particle at the 
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surface position R = (x, y) at time t. This function 
thus generalizes the well known pair distribution function 
,g(R) from statistical mechanics [lH [Hj], since it provides 
information about the interacting particle dynamics. De- 
pending on whether we consider correlations of a particle 
with itself or with other different particles, we distinguish 
between self correlation functions, G S (R, *), and distinct 
correlation functions, Gd(R,*), respectively. With this, 
the full (classical) pair correlation function can be ex- 
pressed as 



G(R,*) = G s (R,*)+G d (R,*). 



(39) 



According to its definition, G S (R, *) is peaked at t = 0, 
and approaches zero as t increases since the particle loses 
correlation with itself. On the other hand, at * = 0, 
Gd(R,t) gives the static pair correlation function (the 
standard pair distribution function), <?(R) = Gd(R, 0), 
while approaches the mean surface number density a of 
diffusing particles as t — > oo. Taking this into account, 
Eq. ([39]) can be expressed as 



G(R, 0) = <J(R) + g(R) 



at * = 0, and as 



G(R,t) S3 a 



(40) 



(41) 



for a homogeneous system with ||R|| — > oo and/or t — > 
oo. At low adparticle concentrations, when interactions 
among adsorbates can be neglected because they are far 
apart from each other, the main contribution to Eq. |39|) 
is G s (particle-particle correlations are negligible and 
Gd ~ 0). On the contrary, for high coverages, it is 
expected that Gd contributes significantly to Eq. ([39]) . 
Within our approach, the interaction among adsorbates 
is described by a particle subjected to a random force (a 
shot noise). Thus, diffusion is described by the G = G s 
function in the interacting single adsorbate approxima- 
tion. At low coverages, the G s function in both approx- 
imations (interacting and noninteracting) has to be the 
same. 

As seen in Sec. m although Gaussian white noise aris- 
ing from the surface is not considered here, an analogous 
Markovian-Langevin also emerges from our model based 
on shot noise. Therefore the same analytical results ob- 
tained elsewhere [6] with a Gaussian white noise can be 
easily extended to our case, where the friction coefficient 
has to be interpreted in terms of the collision frequency 
between adsorbates. For this treatment, the dynamic 
structure factor is better expressed as (22j | 



S(AK,w) 



e -iut ^ e -iAK-R(t) e iAK-R(0)^ df 



J(AK, t) dt, (42) 



where the brackets in the integral denote an ensemble 
average and R(t) is the adparticle trajectory. Here, 



7(AK,f) = ( e -iAK.[R(t)-R(0)]) 



-iAK f* -uk (*') dt'\ 



(43) 



is the intermediate scattering function, which is the space 
Fourier transform of G(R,*) 0]. In Eq. gSJ, u AK 
is the velocity of the adparticle projected onto the di- 
rection of the parallel momentum transfer, with length 
AK = ||AK||. After a second-order cumulant expansion 
in AK in the second equality of Eq. (f43|) , the intermedi- 
ate scattering function reads as 



/(AK,*) 



-AK 2 f*(t-t')C(t')dt' 



(44) 



where C(t) is given by Eq. (I24|) [note that C(t), as given by 
Eq. (f24|) . is independent of the direction of AK]. This is 
the so-called Gaussian approximation 12§ , which is exact 
when the velocity correlations at more than two different 
times are negligible. This allows one to replace the aver- 
age acting over the exponential function in Eq. (|43|) by 
an average acting over its argument, as seen in Eq. (l4"4")l . 

In the case of an almost flat surface the resulting in- 
termediate scattering function is 



/(AK, t) = exp [-x 2 (e~ A * + A* - 1 



with 



X 



)AK 2 /X 2 



(45) 



(46) 



From this relation we can obtain both the mean free 
path, I 

D = T r ( 



Trs/ ( v o)> an d the self-diffusion coefficient, 
(Einstein relation). It can be easily shown 



that the dynamic structure factor derived from Eq. (|45|) 
has a Gaussian shape for short times compared to f = r r 
(or x — > 00 ) an d a Lorentzian shape for long times or 
values of small AK (or x <C 1) @, HH- In this last case, 
Eq. (|45[) becomes a pure decaying exponential function, 



/(AK,*) 



-X 2 At 



(47) 



which is in agreement with the ansatz used in Ref. [3] to 
fit the polarization values obtained experimentally. 

On the other hand, when a harmonic oscillator is con- 
sidered, Eq. (|44l) becomes 



/(AK, t) = exp 



X 2 A 2 

LOLOq 



cos 5 — 



-At/2 



cos(w* — S) 



(48) 

As infers from this equation, the intermediate scatter- 
ing function displays an oscillatory, but exponentially 
damped, behavior around a certain value [although in 
the limit uj — > Eq. (gSJ) approaches Eq. ([IB"]!]. Unlike 
the free-potential case, this means that after relaxation 
some correlation is still present in the system, being the 
limit value 



/(AK,oo) -> e 



~X A cosS/ujujq 



(49) 



Again, this is in agreement with the experimental ob- 
servation. According to Ref. 0] the residual or anoma- 
lous value of the polarization observed seems to be con- 
nected to the blockage of the adsorbate perpendicular 
motion (which is related to translational hopping). On 
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the other hand, according to our model, another source 
for such an effect could be that some Na atoms would 
become trapped inside potential wells, thus behaving as 
damped oscillators. Nonetheless, note that in both cases 
the anomalous behavior for I(t) arises from inhibiting 
the free diffusion of the atoms on the surface. For re- 
alistic corrugated surfaces, one will find a nonseparable 
combination of this behavior (related to the trapped par- 
ticles inside the potential wells) and that described by 
Eq. (|4"5)) (associated to the free-diffusing particles). The 
corresponding scattering law has been given elsewhere 
for models where the velocity autocorrelation func- 
tion can be expressed in terms of a damped anharmonic 
oscillator. 

Taking into account the approximation for the interme- 
diate scattering function, we can now relate in an easy 
manner the velocity power spectrum with the dynamic 
structure factor. As can be readily seen from Eq. 
1231] . we have 



(50) 



where the velocity is recorded along the AK direction. 
The velocity power spectrum or Fourier transform of the 
velocity autocorrelation function is defined as 



Z(w) 



C(t) 



^ luJt dt. 



Thus, by Fourier transforming both sides in Eq. 
obtain 



(51) 



we 



Z{w) 



uj lim 

A if— 



S(AK,cu) 
AK 2 



The interest in this expression relies on the fact that 
it allows us to obtain a relationship between the diffu- 
sion coefficient, D, and the full width at half maximum 
(FWHM) of S by setting u = in Eq. since 



D = -Z(u 
2 v 



0). 



(53) 



For example, for a Lorentzian-shaped S, the FWHM is 
given by T = 2DAK 2 , which is basically the line shape 
that one can observe in an experiment. On the other 
hand, if the line shape is Gaussian, its FWHM is T = 
2v / 2hi2( W 2)i/2 AK. 

Before concluding this section, we would like to give 
a brief account on how to relate in a simple manner the 
coverage 9 and A. In the elementary kinetic theory of 
transport in gases (see, for example, Ref. [HI]), diffusion 
is proportional to the mean free path Z, which is propor- 
tionally inverse to both the density of gas particles and 
the effective area of collision when a hard-sphere model is 
assumed. For 2D collisions, the effective area is replaced 
by an effective length (twice the radius p of the adparti- 
cle) and the gas density by the surface density a. With 
this, the mean free path reads as 



Z = 



1 



2y/2pa' 



(54) 



According to the Chapman-Enskog theory for hard 
spheres, the self-diffusion coefficient can be written as 



D = 



1 

6pa 



(55) 



Now, from the Einstein relation, and taking into account 
that 8 = a 2 a for a square surface lattice of unit cell length 
a, we obtain 



A 



6p9 



knT 



(56) 



Therefore, given a certain surface coverage and temper- 
ature, A can be readily estimated from Eq. (|56p . 



IV. RESULTS AND DISCUSSION 



A. Numerical details 

In order to solve Eq. ([5]) we have used the velocity Ver- 
let algorithm, which is commonly applied when dealing 
with stochastic differential equations 26]. For the aver- 
age calculations shown here a number of 10 000—20 000 
trajectories is sufficient for convergence. The initial con- 
ditions are chosen such that the velocities are distributed 
according to a Maxwell-Boltzmann velocity distribution 
at a temperature T, and the positions follow a uniform 
random distribution along the extension of a single unit 
cell of the potential model used (see below) . Regarding 
(52) the dynamical parameters, we have used A' = 10~ 3 a.u. 
(w 41.3 ps -1 ) and for the mass and radius of the ad- 
particles we have considered those of a Na atom, since 
Na on Cu(001) is a typical system where diffusion has 
been studied and will serve to illustrate our model nu- 
merically. As for the coverage, 6*N a = 1 corresponds to 
one Na atom per Cu(001) surface atom or, equivalently, 
u = 1.53 X 10 19 atom/cm 2 0; a = 2.557 A is the unit 
cell length; and p — 2 A has been used for the atomic 
radius. For example, with these values a collision fric- 
tion A = 5 x 10~ 6 a.u. (« 0.2067 ps^ 1 ) at a surface 
temperature T = 100 K would be related to a coverage 
#Na = 0.059, while the same friction at a temperature 
T = 200 K would be caused by a coverage 6*N a = 0.042. 
One order of magnitude higher for A implies also one or- 
der of magnitude for the coverage at the corresponding 
surface temperatures. 



B. Diffusion in a fiat surface 

First we are going to analyze the case of a flat sur- 
face, which represents fairly well the situation of a low 
corrugated real surface, where the role of the activation 
barrier is negligible. Thus, in Fig. [lja) we observe that, 
effectively, as predicted by the theory, two (time) regimes 
are distinguishable from the mean square displacement 
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FIG. 1: (a) (x' 2 (t)) for two different temperatures and two 
values of the friction coefficient, (b) C(t) for the same tem- 
peratures and friction coefficients as in (a). The temperatures 
used are T = 100 K (thin lines) and T = 200 K (thick lines), 
and the friction coefficients are A = 5 x 10 -6 a.u. (solid lines) 
and A = 5 x 10~ 5 a.u. (dashed lines). 
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FIG. 2: (a) I(t) for two different temperatures and two values 
of the friction coefficient, and AK = 0.88 A -1 , (b) Same as 
(a), but for AK = 1.23 A -1 . The temperatures used are 
T = 100 K (thin lines) and T — 200 K (thick lines), and the 
friction coefficients are A = 5 x 10 -6 a.u. (solid lines) and 
A = 5 x 10~ 5 a.u. (dashed lines). 



(x 2 (t)): parabolic and linear, corresponding to free and 
diffusive motion, respectively. As is apparent, diffusion 
(proportional to the slope of the linear regime) is en- 
hanced by increasing temperature and decreasing the 
friction coefficient. Within the diffusive regime, though 
the presence of many particles tends to inhibit individual 
particle motions, there is still the possibility for particles 
to propagate but at a smaller rate [compare the cuadratic 
increase of (x 2 (t)) during the free motion with its linear 
increase during the diffusive one]. The effect of the dif- 
fusion process can also be seen by looking at the velocity 
autocorrelation function C(t) plotted in Fig. [ljb). As 
seen, the larger A the faster the decay, this result be- 
ing in agreement with Eq. ([28]) as well as its exponential 
decay. The numerical results fit perfectly on this decay, 
but with a slightly lower A than the nominal one because, 
strictly speaking, ours is not a rigorous white shot noise. 
The agreement can also be observed when increasing the 
temperature, leaving A unchanged; the starting value of 
C(t) increases, but its trend, e~ At , remains the same. 

For a flat surface the Gaussian approximation assumed 
in Eq. (|44p for I(t) is exact. This is numerically corrobo- 
rated by looking at the results presented in Fig. ® where 
two different values of AK, 0.88 A" 1 and 1.23 A -1 , are 



analyzed. As is apparent, for both AK values, the ini- 
tial falloff of I(t) displays a Gaussian shape, while for 
longer times it is exponential. This makes that the dy- 
namic structure factor (see Fig. [3]) at small positive (an- 
nihilation events) and negative (creation events) energy 
transfers, i.e., the quasielastic peak region, displays a 
mixed Gaussian-Lorentzian profile describable by using 
the r and incomplete T functions [6|, |25I |. Notice from 
these plots that as A increases the line shape undergoes 
narrowing, unlike what one would expect when having 
collision events. This result for a flat surface, which is 
also predicted analytically, can be explained as follows. 
Noninteracting adsorbates behave like an ideal gas, i.e., 
particles spread out freely, without feeling the action of 
any other particle. This gives rise to dynamic struc- 
ture factors that display Gaussian profiles. However, as 
particle-particle interactions are taken into account there 
is a friction arising from the neighboring adsorbates that 
opposes the free motion, which increases with increas- 
ing A. That is, going back to the results presented in 
Fig- [Ha), the particle reaches the diffusive regime faster 
as A becomes larger. The corresponding line shapes are 
no longer Gaussian functions and display narrowing at 
larger values of A (see Fig. [3]), as analytically predicted 
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FIG. 3: (a) S(ui) for two different temperatures and two values 
of the friction coefficient, and AK = 0.88 A -1 , (b) Same as 
(a), but for AK = 1.23 A -1 . The temperatures used are 
T = 100 K (thin lines) and T = 200 K (thick lines), and 
the friction coefficients are A = 5 x 10~ 6 a.u. (solid lines) 
and A = 5 x 10 _J a.u. (dashed lines). In all cases, S(ui) has 
been normalized to unity in order to better appreciate the 
line shape broadening/narrowing; moreover, they have also 
been smoothed out to eliminate contributions coming from 
long-time fluctuations in I(t). 



by the scattering law [6|, |25( ■ Of course, as temperature 
increases, for a fixed A, one evidently observes the broad- 
ening of the line shape. 

Finally, observe that the results presented here are sim- 
ilar to those that one would obtain with a Brownian-likc 
motion: (i) with A a slow diffusion is observed according 
to the Einstein relation, and (ii) diffusion becomes more 
active as the temperature of the ensemble of adparticles 
increases. However, also note that unlike a Brownian 
motion, diffusion appears here as a consequence of the 
discrete (in time) "kicks" felt by the particles. More- 
over, between two consecutive "kicks" or collisions, they 
move basically without feeling any stochastic force (this 
is similar to the collision model proposed in Ref. [171]). 



C. Diffusion in a separable potential 

When the surface corrugation is relatively strong and 
cannot be neglected, it will induce very important effects 
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FIG. 4: (a) (x 2 (t)} for two different temperatures and two 
values of the friction coefficient, (b) Enlargement at short 
times of part (a), (c) C(t) for the same temperatures and 
friction coefficients as in (a) . The temperatures used are T = 
100 K (thin lines) and T = 200 K (thick lines), and the friction 
coefficients are A = 5 x f0~ 6 a.u. (solid lines) and A = 5 x 
fO -5 a.u. (dashed lines). 



regarding the adsorbate dynamics. To illustrate these 
effects here we consider an adsorbate-surface interaction 
potential model 



V(x, y) — V [2 — cos(27ra;/a) — cos(2ny/a)], 



(57) 



where 2Vb = 33.5 meV is the activation barrier height 
in one direction [x or y). As a function of the temper- 
ature, the behavior is the same as previously observed 
for a flat surface. As can be seen in Fig. SJa) [and in an 
enlargement view in Fig. H^b)] , for a given temperature 
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FIG. 5: (a) I(t) for two different temperatures and two values 
of the friction coefficient, and AK — 0.88 A -1 , (b) Same as 
(a), but for AK = 1.23 A -1 . The temperatures used are: 
T = 100 K (thin lines) and T = 200 K (thick lines); and the 
friction coefficients are: A = 5 x 10~ 6 a.u. (solid lines) and 
A = 5 x 10~ 5 a.u. (dashed lines). 



the presence of the potential also inhibits the diffusion of 
particles at larger values of A since the number of them 
trapped in the potential wells is also larger (according to 
Einstein's law). This result is also reported in Refs. [U, [|| 
when the coverage increases. However, the effect of the 
static interaction potential will give rise to observing the 
opposite behavior than in the case of a flat surface: as 
A increases line shapes become broader, a result which 
is also observed experimentally [3]. This is because now 
both motions, diffusion and vibration, are coupled. Thus, 
with increasing A, the adsorbates can remain localized in- 
side a given surface well for longer times. This allows one 
to observe a higher collisional damping in the T mode, 
which participates actively in the diffusion process. This 
is clearly depicted in Fig. 0|c), where the velocity au- 
tocorrelation function displays an oscillatory behavior, 
which is more damped as A increases. We observe that 
C(t) fits the profile given by Eq. (f3"Tj) but with parame- 
ters different from a harmonic oscillator (see Ref. [6| for 
an anharmonic oscillator model). Only for a motion lo- 
calized mainly at the bottom of the potential well, the 
harmonic oscillator profile will be reproduced. 

The oscillations due to the T mode are not only ob- 
servable in a plot of C(t), but they also manifest in I(t) 
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FIG. 6: (a) S(uj) for two different temperatures and two values 
of the friction coefficient, and AK = 0.88 A -1 , (b) Same as 
(a), but for AK = 1.23 A -1 . The temperatures used are: 
T = 100 K (thin lines) and T = 200 K (thick lines); and 
the friction coefficients are: A = 5 x 10~ 6 a.u. (solid lines) 
and A = 5 x 10~ 5 a.u. (dashed lines). In all cases, S(lo) has 
been normalized to unity in order to better appreciate the 
line shape broadening or narrowing; moreover, they have also 
been smoothed out to eliminate contributions coming from 
long-time fluctuations in I(t). 



(see Fig. [5]) and S(uj) (see Fig. [5]) for two values of AK, 
0.88 A -1 and 1.23 A" 1 , covering the region of the first 
Brillouin zone, [0, 1.3] A -1 . As shown elsewhere @, the 
FWHM of the quasielastic peak is quite similar for both 
AK values and therefore the variation obtained in the 
width is due only to the effect of A. Plotting I(t) we 
can observe a loss of phase when comparing this function 
for different values of A. This behavior continues for all 
oscillations and ends up with I(t) falling faster for the 
case with larger A, thus given rise to broadening as A 
increases. Nonetheless, it is expected that for relatively 
large values of A, one can recover the behavior observed 
in the flat case. 

Finally, in Fig.[7]the peak corresponding to the T mode 
placed around the frequency of oscillation (4 meV in en- 
ergy) is plotted for the same two values of AK as before. 
Not only a shift of the position but also a broadening arc 
clearly seen with A, as also observed experimentally 
For both the quasielastic peak and the T mode peak the 
analytical formulas for the line shapes given elsewhere Q 
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FIG. 7: Same as Fig.[6l but enlarged to show the peaks corre- 
sponding to the T mode. Since this mode is symmetric with 
respect to lu — meV, only the one located to the right of 
the quasielastic peak is shown in every case. 



in the context of Gaussian white noise also fit fairly well 
the numerical results presented for the white shot noise 
model used here. 



V. CONCLUSIONS 



The broadening obtained by this shot noise model 
agrees qualitatively well with the QHAS experimental ob- 
servations [|[ . Note that the broadening described by our 
model is of the same type as that observed in the spec- 
tral lines of gases under high pressure conditions. More- 
over, in our simple model the important issue of how the 
presence of more and more adatoms changes the fields of 
force felt by the remaining adparticles is not addressed. 
As far as we know, this problem has not been treated in 
the literature and certainly merits further investigation. 
It would be very interesting to know how it affects the 
broadening of the quasielastic peak. Nevertheless, at the 
level of the experiments carried out it could happen that 
the statistical limit (a large number of collisions during 
the time scales considered) would not lead to any rele- 
vant feature, since any effect linked to type of interaction 
would be blurred up with time. 

Regarding the validity of our model, since it is purely 
stochastic, in principle there are no other limitations than 
those imposed by realistic physical conditions. For in- 
stance, it is clear that the maximum coverage should be 
#Na = 1- However, as the coverage increases, it is also 
apparent that adsorbate-adsorbate interactions will play 
a more prominent role in the diffusion dynamics, since 
in average the motion of each adsorbate will be slowed 
down and they will feel for longer times the force exerted 
by their neighbors. In such a case, the Markovian ap- 
proximation will break down and memory effects should 
be taken into account. On the other hand, according 
to Ref. 9] there is experimental evidence that for cover- 
ages greater than #Na = 0.05, motion perpendicular to 
the surface is observed, this contributing to the diffusion 
dynamics. Therefore, one might think that a stochastic 
model such as ours could at least work fine up to values 
of the coverage about #Na = 0.2 or greater depending 
on whether QHAS measurements still display effective 
Lorentzian functions for quasielastic line shapes [3[. 



We would like to stress that the same treatment de- 
veloped here for a separable 2D periodic potential can 
be easily generalized to a nonseparable 2D periodic sur- 
face with substrate friction due to a Gaussian white noise 
[l9| . this model describing the process in a more re- 
alistic fashion. In particular, preliminary results with 
the nonseparable 2D periodic surface for a low coverage 
(#Na = 0.028) have rendered a value T = 120 /ieV for the 
quasielastic peak width at AK = 1.26 A -1 , which is in 
agreement with the experimental and theoretical values, 
r cxp = 110 /zeV and r t h = 110 /xeV, respectively, given 
in Ref. 0]. 
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